home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / kcl / akcl / kcl.lha / c / num_sfun.c < prev    next >
C/C++ Source or Header  |  1987-06-04  |  10KB  |  545 lines
  1. /*
  2. (c) Copyright Taiichi Yuasa and Masami Hagiya, 1984.  All rights reserved.
  3. Copying of this file is authorized to users who have executed the true and
  4. proper "License Agreement for Kyoto Common LISP" with SIGLISP.
  5. */
  6.  
  7. #include "include.h"
  8. #include "num_include.h"
  9.  
  10. object imag_unit, minus_imag_unit, imag_two;
  11.  
  12. int
  13. fixnum_expt(x, y)
  14. {
  15.     int z;
  16.  
  17.     z = 1;
  18.     while (y > 0)
  19.         if (y%2 == 0) {
  20.             x *= x;
  21.             y /= 2;
  22.         } else {
  23.             z *= x;
  24.             --y;
  25.         }
  26.     return(z);
  27. }
  28.  
  29. object
  30. number_exp(x)
  31. object x;
  32. {
  33.     double exp();
  34.  
  35.     switch (type_of(x)) {
  36.  
  37.     case t_fixnum:
  38.     case t_bignum:
  39.     case t_ratio:
  40.         return(make_shortfloat((shortfloat)exp(number_to_double(x))));
  41.  
  42.     case t_shortfloat:
  43.         return(make_shortfloat((shortfloat)exp((double)(sf(x)))));
  44.  
  45.     case t_longfloat:
  46.         return(make_longfloat(exp(lf(x))));
  47.  
  48.     case t_complex:
  49.     {
  50.         object y, y1;
  51.         object number_sin(), number_cos();
  52.             vs_mark;
  53.     
  54.         y = x->cmp.cmp_imag;
  55.         x = x->cmp.cmp_real;
  56.         x = number_exp(x);
  57.         vs_push(x);
  58.         y1 = number_cos(y);
  59.         vs_push(y1);
  60.         y = number_sin(y);
  61.         vs_push(y);
  62.         y = make_complex(y1, y);
  63.         vs_push(y);
  64.         x = number_times(x, y);
  65.         vs_reset;
  66.         return(x);
  67.     }
  68.  
  69.     default:
  70.         FEwrong_type_argument(Snumber, x);
  71.     }
  72. }
  73.  
  74. object
  75. number_expt(x, y)
  76. object x, y;
  77. {
  78.     enum type tx, ty;
  79.     object z, number_nlog();
  80.     vs_mark;
  81.  
  82.     tx = type_of(x);
  83.     ty = type_of(y);
  84.     if (ty == t_fixnum && fix(y) == 0)
  85.         switch (tx) {
  86.         case t_fixnum:  case t_bignum:  case t_ratio:
  87.             return(small_fixnum(1));
  88.  
  89.         case t_shortfloat:
  90.             return(make_shortfloat((shortfloat)1.0));
  91.  
  92.         case t_longfloat:
  93.             return(make_longfloat(1.0));
  94.  
  95.         case t_complex:
  96.             z = number_expt(x->cmp.cmp_real, y);
  97.             vs_push(z);
  98.             z = make_complex(z, small_fixnum(0));
  99.             vs_reset;
  100.             return(z);
  101.  
  102.         default:
  103.             FEwrong_type_argument(Snumber, x);
  104.         }
  105.     if (number_zerop(x)) {
  106.         if (!number_plusp(ty==t_complex?y->cmp.cmp_real:y))
  107.             FEerror("Cannot raise zero to the power ~S.", 1, y);
  108.         return(number_times(x, y));
  109.     }
  110.     if (ty == t_fixnum || ty == t_bignum) {
  111.         if (number_minusp(y)) {
  112.             z = number_negate(y);
  113.             vs_push(z);
  114.             z = number_expt(x, z);
  115.             vs_push(z);
  116.             z = number_divide(small_fixnum(1), z);
  117.             vs_reset;
  118.             return(z);
  119.         }
  120.         z = small_fixnum(1);
  121.         vs_push(z);
  122.         vs_push(Cnil);
  123.         vs_push(Cnil);
  124.         while (number_plusp(y))
  125.             if (number_evenp(y)) {
  126.                 x = number_times(x, x);
  127.                 vs_top[-1] = x;
  128.                 y = integer_divide1(y, small_fixnum(2));
  129.                 vs_top[-2] = y;
  130.             } else {
  131.                 z = number_times(z, x);
  132.                 vs_top[-3] = z;
  133.                 y = number_minus(y, small_fixnum(1));
  134.                 vs_top[-2] = y;
  135.             }
  136.         vs_reset;
  137.         return(z);
  138.     }
  139.     z = number_nlog(x);
  140.     vs_push(z);
  141.     z = number_times(z, y);
  142.     vs_push(z);
  143.     z = number_exp(z);
  144.     vs_reset;
  145.     return(z);
  146. }
  147.  
  148. object
  149. number_nlog(x)
  150. object x;
  151. {
  152.     double log();
  153.     object r, i, a, p, number_sqrt(), number_atan2();
  154.     vs_mark;
  155.  
  156.     if (type_of(x) == t_complex) {
  157.         r = x->cmp.cmp_real;
  158.         i = x->cmp.cmp_imag;
  159.         goto COMPLEX;
  160.     }
  161.     if (number_zerop(x))
  162.         FEerror("Zero is the logarithmic singularity.", 0);
  163.     if (number_minusp(x)) {
  164.         r = x;
  165.         i = small_fixnum(0);
  166.         goto COMPLEX;
  167.     }
  168.     switch (type_of(x)) {
  169.     case t_fixnum:
  170.     case t_bignum:
  171.     case t_ratio:
  172.         return(make_shortfloat((shortfloat)log(number_to_double(x))));
  173.  
  174.     case t_shortfloat:
  175.         return(make_shortfloat((shortfloat)log((double)(sf(x)))));
  176.  
  177.     case t_longfloat:
  178.         return(make_longfloat(log(lf(x))));
  179.  
  180.     default:
  181.         FEwrong_type_argument(Snumber, x);
  182.     }
  183.  
  184. COMPLEX:
  185.     a = number_times(r, r);
  186.     vs_push(a);
  187.     p = number_times(i, i);
  188.     vs_push(p);
  189.     a = number_plus(a, p);
  190.     vs_push(a);
  191.     a = number_nlog(a);
  192.     vs_push(a);
  193.     a = number_divide(a, small_fixnum(2));
  194.     vs_push(a);
  195.     p = number_atan2(i, r);
  196.     vs_push(p);
  197.     x = make_complex(a, p);
  198.     vs_reset;
  199.     return(x);
  200. }
  201.  
  202. object
  203. number_log(x, y)
  204. object x, y;
  205. {
  206.     object z;
  207.     vs_mark;
  208.  
  209.     if (number_zerop(y))
  210.         FEerror("Zero is the logarithmic singularity.", 0);
  211.     if (number_zerop(x))
  212.         return(number_times(x, y));
  213.     x = number_nlog(x);
  214.     vs_push(x);
  215.     y = number_nlog(y);
  216.     vs_push(y);
  217.     z = number_divide(y, x);
  218.     vs_reset;
  219.     return(z);
  220. }
  221.  
  222. object
  223. number_sqrt(x)
  224. object x;
  225. {
  226.     object z;
  227.     double sqrt();
  228.     vs_mark;
  229.  
  230.     if (type_of(x) == t_complex)
  231.         goto COMPLEX;
  232.     if (number_minusp(x))
  233.         goto COMPLEX;
  234.     switch (type_of(x)) {
  235.     case t_fixnum:
  236.     case t_bignum:
  237.     case t_ratio:
  238.         return(make_shortfloat(
  239.             (shortfloat)sqrt(number_to_double(x))));
  240.  
  241.     case t_shortfloat:
  242.         return(make_shortfloat((shortfloat)sqrt((double)(sf(x)))));
  243.  
  244.     case t_longfloat:
  245.         return(make_longfloat(sqrt(lf(x))));
  246.  
  247.     default:
  248.         FEwrong_type_argument(Snumber, x);
  249.     }
  250.  
  251. COMPLEX:
  252.     z = make_ratio(small_fixnum(1), small_fixnum(2));
  253.     vs_push(z);
  254.     z = number_expt(x, z);
  255.     vs_reset;
  256.     return(z);
  257. }
  258.  
  259. object
  260. number_atan2(y, x)
  261. object y, x;
  262. {
  263.     object z;
  264.     double atan(), dy, dx, dz;
  265.  
  266.     dy = number_to_double(y);
  267.     dx = number_to_double(x);
  268.     if (dx > 0.0)
  269.         if (dy > 0.0)
  270.             dz = atan(dy / dx);
  271.         else if (dy == 0.0)
  272.             dz = 0.0;
  273.         else
  274.             dz = -atan(-dy / dx);
  275.     else if (dx == 0.0)
  276.         if (dy > 0.0)
  277.             dz = PI / 2.0;
  278.         else if (dy == 0.0)
  279.             FEerror("Logarithmic singularity.", 0);
  280.         else
  281.             dz = -PI / 2.0;
  282.     else
  283.         if (dy > 0.0)
  284.             dz = PI - atan(dy / -dx);
  285.         else if (dy == 0.0)
  286.             dz = PI;
  287.         else
  288.             dz = -PI + atan(-dy / -dx);
  289.     if (type_of(x) == t_longfloat || type_of(y) == t_longfloat)
  290.         z = make_longfloat(dz);
  291.     else
  292.         z = make_shortfloat((shortfloat)dz);
  293.     return(z);
  294. }
  295.  
  296. object
  297. number_atan(y)
  298. object y;
  299. {
  300.     object z, z1;
  301.         vs_mark;
  302.  
  303.     if (type_of(y) == t_complex) {
  304.         z = number_times(imag_unit, y);
  305.         vs_push(z);
  306.         z = one_plus(z);
  307.         vs_push(z);
  308.         z1 = number_times(y, y);
  309.         vs_push(z1);
  310.         z1 = one_plus(z1);
  311.         vs_push(z1);
  312.         z1 = number_sqrt(z1);
  313.         vs_push(z1);
  314.         z = number_divide(z, z1);
  315.         vs_push(z);
  316.         z = number_nlog(z);
  317.         vs_push(z);
  318.         z = number_times(minus_imag_unit, z);
  319.         vs_reset;
  320.         return(z);
  321.     }
  322.     return(number_atan2(y, small_fixnum(1)));
  323. }
  324.  
  325. object
  326. number_sin(x)
  327. object x;
  328. {
  329.     double sin();
  330.  
  331.     switch (type_of(x)) {
  332.  
  333.     case t_fixnum:
  334.     case t_bignum:
  335.     case t_ratio:
  336.         return(make_shortfloat((shortfloat)sin(number_to_double(x))));
  337.  
  338.     case t_shortfloat:
  339.         return(make_shortfloat((shortfloat)sin((double)(sf(x)))));
  340.  
  341.     case t_longfloat:
  342.         return(make_longfloat(sin(lf(x))));
  343.  
  344.     case t_complex:
  345.     {
  346.         object    r;
  347.         object    x0, x1, x2;
  348.         vs_mark;
  349.  
  350.         x0 = number_times(imag_unit, x);
  351.         vs_push(x0);
  352.         x0 = number_exp(x0);
  353.         vs_push(x0);
  354.         x1 = number_times(minus_imag_unit, x);
  355.         vs_push(x1);
  356.         x1 = number_exp(x1);
  357.         vs_push(x1);
  358.         x2 = number_minus(x0, x1);
  359.         vs_push(x2);
  360.         r = number_divide(x2, imag_two);
  361.  
  362.         vs_reset;
  363.         return(r);
  364.     }
  365.  
  366.     default:
  367.         FEwrong_type_argument(Snumber, x);
  368.  
  369.     }
  370. }
  371.  
  372. object
  373. number_cos(x)
  374. object x;
  375. {
  376.     double cos();
  377.  
  378.     switch (type_of(x)) {
  379.  
  380.     case t_fixnum:
  381.     case t_bignum:
  382.     case t_ratio:
  383.         return(make_shortfloat((shortfloat)cos(number_to_double(x))));
  384.  
  385.     case t_shortfloat:
  386.         return(make_shortfloat((shortfloat)cos((double)(sf(x)))));
  387.  
  388.     case t_longfloat:
  389.         return(make_longfloat(cos(lf(x))));
  390.  
  391.     case t_complex:
  392.     {
  393.         object r;
  394.         object x0, x1, x2;
  395.         vs_mark;
  396.  
  397.         x0 = number_times(imag_unit, x);
  398.         vs_push(x0);
  399.         x0 = number_exp(x0);
  400.         vs_push(x0);
  401.         x1 = number_times(minus_imag_unit, x);
  402.         vs_push(x1);
  403.         x1 = number_exp(x1);
  404.         vs_push(x1);
  405.         x2 = number_plus(x0, x1);
  406.         vs_push(x2);
  407.         r = number_divide(x2, small_fixnum(2));
  408.  
  409.         vs_reset;
  410.         return(r);
  411.     }
  412.  
  413.     default:
  414.         FEwrong_type_argument(Snumber, x);
  415.  
  416.     }
  417. }
  418.  
  419. object
  420. number_tan(x)
  421. object x;
  422. {
  423.     object r, s, c;
  424.     vs_mark;
  425.  
  426.     s = number_sin(x);
  427.     vs_push(s);
  428.     c = number_cos(x);
  429.     vs_push(c);
  430.     if (number_zerop(c) == TRUE)
  431.         FEerror("Cannot compute the tangent of ~S.", 1, x);
  432.     r = number_divide(s, c);
  433.     vs_reset;
  434.     return(r);
  435. }
  436.  
  437. Lexp()
  438. {
  439.     check_arg(1);
  440.     check_type_number(&vs_base[0]);
  441.     vs_base[0] = number_exp(vs_base[0]);
  442. }
  443.  
  444. Lexpt()
  445. {
  446.     check_arg(2);
  447.     check_type_number(&vs_base[0]);
  448.     check_type_number(&vs_base[1]);
  449.     vs_base[0] = number_expt(vs_base[0], vs_base[1]);
  450.     vs_pop;
  451. }
  452.  
  453. Llog()
  454. {
  455.     int narg;
  456.     
  457.     narg = vs_top - vs_base;
  458.     if (narg < 1)
  459.         too_few_arguments();
  460.     else if (narg == 1) {
  461.         check_type_number(&vs_base[0]);
  462.         vs_base[0] = number_nlog(vs_base[0]);
  463.     } else if (narg == 2) {
  464.         check_type_number(&vs_base[0]);
  465.         check_type_number(&vs_base[1]);
  466.         vs_base[0] = number_log(vs_base[1], vs_base[0]);
  467.         vs_pop;
  468.     } else
  469.         too_many_arguments();
  470. }
  471.  
  472. Lsqrt()
  473. {
  474.     check_arg(1);
  475.     check_type_number(&vs_base[0]);
  476.     vs_base[0] = number_sqrt(vs_base[0]);
  477. }
  478.  
  479. Lsin()
  480. {
  481.     check_arg(1);
  482.     check_type_number(&vs_base[0]);
  483.     vs_base[0] = number_sin(vs_base[0]);
  484. }
  485.  
  486. Lcos()
  487. {
  488.     check_arg(1);
  489.     check_type_number(&vs_base[0]);
  490.     vs_base[0] = number_cos(vs_base[0]);
  491. }
  492.  
  493. Ltan()
  494. {
  495.     check_arg(1);
  496.     check_type_number(&vs_base[0]);
  497.     vs_base[0] = number_tan(vs_base[0]);
  498. }
  499.  
  500. Latan()
  501. {
  502.     int narg;
  503.  
  504.     narg = vs_top - vs_base;
  505.     if (narg < 1)
  506.         too_few_arguments();
  507.     if (narg == 1) {
  508.         check_type_number(&vs_base[0]);
  509.         vs_base[0] = number_atan(vs_base[0]);
  510.     } else if (narg == 2) {
  511.         check_type_or_rational_float(&vs_base[0]);
  512.         check_type_or_rational_float(&vs_base[1]);
  513.         vs_base[0] = number_atan2(vs_base[0], vs_base[1]);
  514.         vs_pop;
  515.     } else
  516.         too_many_arguments();
  517. }
  518.  
  519. init_num_sfun()
  520. {
  521.     imag_unit
  522.     = make_complex(make_shortfloat((shortfloat)0.0),
  523.                make_shortfloat((shortfloat)1.0));
  524.     enter_mark_origin(&imag_unit);
  525.     minus_imag_unit
  526.     = make_complex(make_shortfloat((shortfloat)0.0),
  527.                make_shortfloat((shortfloat)-1.0));
  528.     enter_mark_origin(&minus_imag_unit);
  529.     imag_two
  530.     = make_complex(make_shortfloat((shortfloat)0.0),
  531.                make_shortfloat((shortfloat)2.0));
  532.     enter_mark_origin(&imag_two);
  533.  
  534.     make_constant("PI", make_longfloat(PI));
  535.  
  536.     make_function("EXP", Lexp);
  537.     make_function("EXPT", Lexpt);
  538.     make_function("LOG", Llog);
  539.     make_function("SQRT", Lsqrt);
  540.     make_function("SIN", Lsin);
  541.     make_function("COS", Lcos);
  542.     make_function("TAN", Ltan);
  543.     make_function("ATAN", Latan);
  544. }
  545.